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ABSTRACT 

Most of the super massive black hole mass (M.) estimates based on stellar kinematics use 
the assumption that galaxies are axisymmetric oblate spheroids or spherical. Here we use 
fully general triaxial orbit-based models to explore the effect of relaxing the axisymmetric 
assumption on the previously studied galaxies M32 and NGC 3379. We find that M32 can only 
be modeled accurately using an axisymmetric shape viewed nearly edge-on and our black 
hole mass estimate is identical to previous studies. When the observed 5° kinematical twist 
is included in our model of NGC 3379, the best shape is mildly triaxial and we find that our 
best-fitting black hole mass estimate doubles with respect to the axisymmetric model. This 
particular black hole mass estimate is still within the errors of that of the axisymmetric model 
and consistent with the M.-a relationship. However, this effect may have a pronounced impact 
on black hole demography, since roughly a third of the most massive galaxies are strongly 
triaxial. 

Key words: black hole physics, galaxies: elliptical and lenticular, cD - galaxies: kinematics 
and dynamics - galaxies: individual: NGC 3379, M32 galaxies: structure, galaxies: nuclei 



1 INTRODUCTION 

The masses of super massive black holes in the centers of galaxies 
are known to correlate with several properties of the host galaxy. The 
most well known correlation is with the stellar velocity dispersion 
of the galaxy (M.-a, e.g., Tremaine et al. 2002). The black hole is 
thought to play an important role in the evolution of its host (e.g. 
AGN feedback), as the properties of galaxies are tightly linked to M.. 
Therefore it is important to be able to measure the M. accurately and 
understand whether the scatter in the relations is due to measurement 
error, or if it is intrinsic. Most of the M. estimates upon which these 
relationships are based were derived using dynamical (edge-on) 
axisymmetric (or spherical) dynamical models (See Ferrarese & 
Ford 2005 for a review). 

It has long been known from photometry that some elliptical 
galaxies are triaxial (e.g. Kormendy & Bender 1996), i.e. have 
intrinsic shapes with three distinct axes (Binney 1976, 1978), and 
more recently Emsellem et al. (2007) and Cappellari et al. (2007) 
have shown using stellar kinematics that around a third of the most 
massive ellipticals are at least mildly triaxial. It is thus relevant 
to rederive the M. in these galaxies with our triaxial instead of 
axisymmetric orbit super-position models. 

Thomas et al. (2007) have modelled mock triaxial merger rem- 
nants using axisymmetric geometry and found a correlation between 
viewing angle and the recovered total galaxy mass. Additionally, 
they found that the luminous mass to light ratio (M/L) was underesti- 
mated by up to 20% confirming that triaxial modelling is preferred. 

In this paper we explore the black hole recovery with the triax- 



ial machinery from van den Bosch et al. (2008) using galaxies that 
have previously been modeled with axisymmetric codes to verify 
the triaxial models. We do this with M32 and NGC 3379 which have 
their black hole mass determined using STIS, as well as SAURON and 
OASIS Integral Field Unit (IFU) data. Both galaxies have state-of- 
the-art kinematics available over a large spatial range and are inside 
the sphere of influence of the black hole. 

We describe our modeling technique and uncertainties in §2 and 
then derive our black hole estimates on galaxies M32 and NGC 3379 
in §3. We briefly address the reliability in §4, and we end with 
discussion and conclusions in §5. 



2 TRIAXIAL SCHWARZSCHILD MODELING 

In this paper we use the triaxial Schwarzschild (1979) orbit superpo- 
sition technique as it is described in van den Bosch et al. (2008), of 
which we give a brief summary here. It is a powerful tool to construct 
realistic dynamical models. It allows for an arbitrary triaxial gravita- 
tional potential (with possible contributions from dark components) 
in which the equations of motion are integrated numerically for a 
representative library of orbits. Then the superposition of orbits is 
determined for which the combined density and velocity moments 
best fit the observed surface brightness and kinematics using least 
squares. By marginalising over parameter space, Schwarzschild's 
method not only provides the viewing direction and the M/L with 
the dark matter contribution, but also allows the investigation of the 
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intrinsic dynamical structures as well as the distribution function 
through the orbital mass weights (cf. Vandervoort 1984). These mod- 
els have complete freedom: specifically no form of (an-)isotropy is 
implied, within the limits of the observed photometry and stellar 
kinematics. 

Krajnovic et al. (2005) showed that the (31) Schwarzschild 
method can recover the phase-space distribution function in the 
two-integral axisymmetric case. Additionally, van de Ven, de Zeeuw 
& van den Bosch (2008) confirmed that the same is true for triax- 
ial Schwarzschild models and extended this to show that we can 
recover the internal dynamics and three-integral distribution func- 
tion of triaxial early-type galaxies given their viewing angles. In 
van den Bosch & van de Ven (2009) we have done tests on several 
three-integral Abel models, that simulate real early-type galaxies, 
and found we are also able to recover the viewing angles of triaxial 
early-type galaxies, as long as their kinematics show clearly defined 
gradients. These final tests allowed us to firmly establish robust- 
ness of the shape recovery and phase space distribution of triaxial 
Schwarzschild modeling. 

To construct a (luminous) mass model we assume that the three- 
dimensional mass distribution can be parameterized with multiple 
coaxial Gaussians (Monnet et al. 1992; Bendinelli & Parmeggiani 
1995; Emsellem et al. 1994; Cappellari 2002). The mass distribution 
used in the models is constructed by fitting two-dimensional Gaus- 
sians to the broad band photometry of the galaxy. These Gaussians 
can be deprojected onto a coaxial triaxial shape by choosing three 
viewing angles; (■&, (p) and the apparent misalignment (\|/), which are 
used to define the direction from which the galaxy is seen. The shape 
of each deprojected Gaussian depends on the viewing angles, its 
observed flattening and its isophotal twist. The shape of each three- 
dimensional Gaussian is characterized by pj = bj/cij, qj = cj/cij 
and uj = a'j/aj, where aj, bj and Cj are the long, intermediate and 
short axis lengths of each individual Gauss j, and a'j is the length 
of longest axis of the Gauss as observed on the sky. To do a full 
search of all possible mass models we first sample uniformly over 
a (separate) set of axis ratios, p, q and w, which translates to a set 
of viewing angles. These viewing angles are then used to construct 
all the mass models. The light model is converted to a mass by 
assuming a constant M/L ratio (but see van den Bosch et al. 2006; 
van de Ven et al. 2006). 

The black hole mass estimates determined using Schwarzschild 
modeling have not been without debate: Valluri et al. (2004) reported 
the existence of a % 2 plateau prohibiting a M. assessment. This could 
be avoided by using some form of regularisation and enough orbits 
in the models. The debate was settled by Magorrian (2006) who 
showed that the results from standard superposition methods are 
accurate, if observational errors are taken into account. Currently, 
McDermid is leading a joint effort in a comparative study of the 
M. recovery using three independent axisymmetric codes. Using 
generalised cross validation, he finds that regularisation should be 
used to avoid the known issue of these models over-fitting the data, 
and thus yield reliable error estimates on M. . While these results are 
obtained with axisymmetric modeling, it is likely that they will hold 
in the triaxial case too, as the triaxial method is conceptually very 
similar. In this paper we repeated all modeling with and without 
regularisation and found no difference in the recovered black hole 
masses. 

To determine the uncertainty on the derived shape we shall 
use the x 2 based confidence interval that was established in van 
den Bosch & van de Ven (2009). The intervals are based upon 
the expected standard deviation of the A% 2 (=^/2N i, s , where N„t, s 
is the number of kinematical observations used to constrain the 



model). As the M. determination is only influenced by a few of the 
innermost kinematical observations, we shall use the standard formal 
la and 3o results based upon a % 2 distribution with two degrees of 
freedom. This is based on the assumption that the determined shape 
is independent of the M.. It is not obvious that this assumption holds, 
but as we shall see in the rest of this paper, it works and we discuss 
the validity of this assumption in §4. 



3 BLACK HOLE ESTIMATES USING TRIAXIAL 
MODELS 

Here we describe the results of our dynamical modelling (as 
described in §2) of M32 and NGC 3379. We compare the near- 
axisymmetric M. estimates with literature values and explore the 
effect of the possible triaxial shape upon the black hole mass esti- 
mate and orbital structure. 

3.1 M32 

To be able to compare our results directly to other studies it was 
important that we used galaxies that have a published M. obtained 
using axisymmetric Schwarzschild models. 

Therefore we chose the nearby compact fast rotator E3 galaxy 
M32 (Kormendy 1985), as it has been investigated by several in- 
dependent authors (e.g. van der Marel et al. 1998; Verolme et al. 
2002; Kormendy 2004; Valluri et al. 2004). It is consistent with ax- 
isymmetry, as it shows regular rotation and has almost no isophotal 
twist (< 3°, Peletier 1993; Lauer et al. 1998). For our modeling 
of M32 we assumed a distance of 0.79 Mpc and used the surface 
brightness distribution, and the wide field SAURON data from Cap- 
pellari et al. (2007) and the STIS data from Joseph et al. (2001), 
that probes well within the sphere of influence of the black hole. 
The total number of kinematical observations used to constrain the 
models is 58 from STIS and 964 from the SAURON observations, 
measured up to the Gauss-Hermite moment (van der Marel & Franx 
1993; Gerhard 1993; Rix et al. 1997) h 4 . The SAURON kinematics 
were point-symmetrized as described in Appendix A. Because M32 
has a dispersion (~90 km s _1 ) below the instrumental dispersion of 
SAURON (120 kms -1 ) the moments /13 and higher are hard to mea- 
sure, and thus have large associated errors (Cappellari & Emsellem 
2004; Emsellem et al. 2004). However, they are still included in the 
models to ensure that the (otherwise unconstrained) reconstructed 
LOSVDs do not deviate too far from a Gaussian-like shape. 

We did not use the high spatial resolution (HR) SAURON obser- 
vations from Verolme et al. (2002), as the kinematic extraction of 
the low resolution data by Cappellari et al. (2007) is superior, due to 
a much better extraction method and significantly improved stellar 
templates. Since we have high resolution data from STIS, we chose 
not to re-reduce the HR data. 

3.1.1 Triaxial models of M32 

Before we can get an M. estimate we need to explore the shape 
of M32, because it is too expensive to make models and explore 
parameter space for the full range in shape and M., as it would 
increase the computation time by a factor 10. Therefore, we first 
fix the M. at 2.6 x 1O 6 M , while varying the shape and the M/L. 
By doing this we assume that the derived shape does not depend 
on the specific M. we use. We show in §3.1.2 and §4 that this is a 
reasonable assumption. 

The result of this triaxial modeling is shown in Fig. 1. The 
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Figure 1. Contour maps of the confidence intervals for the parameters of 
the mass model of the models of M32. From blue to red the color denote 
high to low confidence. Areas for which the surface brightness cannot be 
deprojected are white. The contours denote 1, 3 (thick line) sigma confidence 
levels at A% 2 = 63 and 189. The six upper-right panels show the intrinsic 
shape parameters (p,q.u) and the M/L; the three lower left panels show the 
viewing angles (■&, <p, . The combination of r> and (p is shown in a Lambert 
equal area projection (The half circle), seen down the north pole (z-axis). 
In this panel the x, y and z symbols give the location of views down those 
axis. The best-fitting models are very round and near oblate. The preferred 
viewing angles are near the y-axis with nearly no intrinsic misalignment. See 
§2 and §3.1.1 



best-fitting shape is nearly as round as the observed flattening allows 
and is consistent with an oblate axisymmetric spheroid with axis 
ratios (p,q) = (0.95±0.05,0.76tj]^). This is fully consistent with 
expectations given the aligned bi-symmetric rotation in the observed 
velocity field. Curiously though, Verolme et al. (2002) did find a 
best-fitting inclination 70° ± 10° (equivalent to p=l, g=0.73±0.03) 
using axisymmetric modeling. Our results do agree, but our errors 
are different, due to our conservative error bars. Our 'axisymmetric' 
models are still slightly triaxial (p > 0.99) and thus do allow for 
additional freedom from the triaxial orbital families. The box orbits 
are important as our (near) axisymmetric models contain 8% box 
orbits within one R e , 

Within the subset of axisymmetric models, the inclination is not 
well constrained at iJ = i =90° Essentially the full inclination 
range (given the observed flattening) is allowed at the 3a level. 
This result is interesting as it agrees with the observation made in 
Krajnovic et al. (2005), and in van den Bosch & van de Ven (2009), 
that the axisymmetric models cannot constrain the inclination well. 
However, as is also shown in van den Bosch & van de Ven (2009), 
we can marginalize over the allowed shapes to gain insight into the 
physical parameters - like M/L and anisotropy - which tells us more 
about the galaxy than the inclination. The M/L is 1.4 ±0.2 Mq/Lqj, 
identical to the value from the models of Cappellari et al. (2007) 1 



1 Cappellari et al. (2006) finds 1.2 M /L 0i /. 
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Figure 2. M. (top) and the M/L (bottom) confidence levels as a function of 
shape for M32. The small red crosses indicate where models were computed 
and the big black cross indicates the best-fit model. From far left to far right 
the models change from prolate (p = 0.8, q = 0.76), triaxial (p = 0.9, q = 
0.76) and oblate (p = 1.0, q = 0.46). The vertical line denotes the transition 
from triaxial to oblate. The models on the vertical line are the roundest 
possible models. The contours show increasing level of confidence. The inner 
and outer thicrk contour indicate a A% 2 of 14.1 and 63, which correspond 
to the 99.9% and 96% confidence on the M. and shape, respectively. The 
models with M.= are relocated to 8x 10 5 Mq, to place them within the 
plot. 

and Verolme et al. (2002) (when corrected for the different assumed 
distance and for galactic reddening). 

3.1.2 M32 M. estimate 

To illustrate the effect of triaxiality, we show what happens when we 
sample different shapes while varying the M. and the M/L. We vary 
the shape from maximally prolate to maximally oblate (6 models) 
and at different axisymmetric shapes (8 models). These shapes 
continuously follow the q=Q.16 and then the p=\ line in the upper 
leftmost plot (p vs. q) in Fig. 1. This also fully encompasses the 
uncertainty in estimated shape. 

In Fig. 2 we show the confidence levels (A% 2 ) of the M/L and 
M. as a function of these shapes. As expected the contours show 
that the roundest models are the best-fit (g>0.66, p>0.88). After 
marginalizing over the shapes, our M. mass estimate of (2. 4 ± 1.0) x 
10 6 M is fully consistent with previous results of (2 — 4) x 10 6 M e 
from Joseph et al. (2001), (2.4 ± 0.7) x 10 6 M from van der Marel 
et al. (1998) and (2.5 ±0.5) x 10 6 M Verolme et al. (2002) and 
2.6 x 10 6 M Q from M.-o (Tremaine et al. 2002). 

The two most nearly prolate models in this test are a signifi- 
cantly worse fit than the best-fit model. To illustrate this we show 
the kinematics of the observed stellar kinematics and four models, 
varying from oblate to prolate, in Fig. 3. While the oblate and round 
models reproduce the observations, the prolate models are unable to 
do so, and we can therefore conclude that M32 does not have prolate 
shape. 
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Figure 3. Observed stellar kinematics of M32 (left) and the models that correspond to Fig. 2, see labels. Top row is the stellar mean velocity and bottom is the 
stellar velocity dispersion. From left to right: SAURON observations, models from oblate to prolate; (p, u) = (1.00,0.48), (0.97,0.76), (0.84,0.76), (0.77,0.76). 
Contours show representative isophotes. All but the prolate models reproduce the observations well. 



Interestingly the best-fitting M/L and M. do not change sig- 
nificantly for the individual triaxial models (See Fig. 2), indicating 
that our initial assumption, i.e. that the recovered shape is not de- 
pendent on the initial M., is probably reasonable in this case. To 
ensure that our models do not depend on the number of orbits, we 
doubled the number of orbits, and repeated the above test and found 
no significant change on the best-fit error and the formal error bars. 

Overall, this shows that it is possible to recover M. with our 
triaxial method. Furthermore, we showed that, in this case, the 
recovered M. does not significantly depend on the intrinsic shape, 
and that the shape of M32 is very close to oblate. 

3.2 NGC3379 

As the second test galaxy we chose the fast rotator NGC 3379, which 
also has its black hole mass measured with orbit-based models 
(Gebhardt et al. 2000, Shapiro et al. 2006, hereafter S06, Douglas 
et al. 2007), has a decoupled nuclear gas ring (e.g. Statler 2001) and 
an M. estimate from the gas kinematics (S06). The velocity maps 
show regular rotation that is consistent with oblate axisymmetry, and 
probe to well within the sphere of influence. Detailed axisymmetric 
orbit-based and gas disc models are shown in S06. However this 
galaxy shows mild evidence for triaxiality: It has a small isophotal 
twist (Capaccioli et al. 1987), an 5 ± 3° kinematical misalignment 
(S06, Statler et al. 1999) and even shows hints of kinematical twist 
(Krajnovic et al. 2008). Capaccioli et al. (1991) suggested it might 
be a triaxial SO seen face-on, and Statler (2001) has argued that 
this galaxy is mildly triaxial. Krajnovic et al. (2008) showed that 
NGC 3379 is consistent with being like all the other fast-rotators in 
the SAURON sample: nearly axisymmetric and with a disk, but seen 
almost face-on (based on their V/o diagram). Their interpretation 
would mean that the misalignment is due to a non-perfect circularity 
of the disk. 

All this makes NGC 3379 an ideal test case for the black hole re- 
covery of the triaxial Schwarzschild machinery. We use the surface 



brightness distribution (based upon WFPC2 and ground-based imag- 
ing), SAURON and OASIS stellar kinematics, and distance (10.28 
Mpc) from S06 for our modeling. 



3.2.1 NGC 3379 in the axisymmetric limit 

NGC 3379 is a very round galaxy, so it is difficult to establish the 
photometric PA with high precision, and thus the exact amount of 
misalignment is uncertain (5 ± 3°), and is in principle consistent 
with no misalignment (0°) within 2a. As such this galaxy could be 
an axisymmetric oblate spheroid. Exploiting this uncertainty, S06 
aligned the stellar kinematics with the photometry by rotating the 
velocity field by 5° and bi-symmetrizing the kinematics. This adjust- 
ment to the observed kinematics ensures that they are completely 
compatible with the axisymmetric modelling and these 'corrected' 
observations were used for their M. determinations. To check our 
triaxial modeling machinery, we duplicate those models from S06 
to see if we recover the same M. given the same input parameters. 
This test is important as our triaxial machinery is not capable of 
making a perfectly axisymmetric model, so even our axisymmetric 
model would be minutely triaxial (p > 0.99) and can still benefit 
from the additional orbital families, especially the strongly radial 
box orbits, which do not exist in a pure axisymmetric potential. 

We reconstructed the parameter space in M/L and M. from 
S06. The results, presented in Fig. 4, are identical to the estimate 
produced with the axisymmetric code, giving M.=(1.4±0.9) x 10 8 
M0 and an M/L of (3.1 ±0.2) Mq/Lqj and is also similar to the 
results from Gebhardt et al. (2000). The fit to the kinematics is 
not shown, as it is essentially identical to the figures in S06. To 
be completely confident we also tested with double the number of 
orbits and this does not affect the recovered black hole mass and its 
confidence interval. 
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Figure 4. The confidence levels in the M. and the M/L for an axisymmetric 
model of NGC 3379 made with the triaxial software. See §3.2.1. The thick 
contour denotes Ax 2 =14.1, corresponding to a 99.9% confidence level. The 
location of the minimum and the shape of the confidence contours match the 
results from S06 (see the rightmost panel in their Fig. 8), demonstrating that 
both codes provide identical results when given the same observables and 
shape. 

3.2.2 Triaxial models of NGC 3379 

To study the effect of kinematical misalignment on the M. esti- 
mate we modelled NGC 3379 using our triaxial machinery and with 
the kinematical misalignment intact and point-symmetrized (See 
Appendix A) kinematics. We first searched through the shape dis- 
tribution, while keeping M. fixed at 1.4 x 10 s M . The best-fitting 
kinematics is shown in Fig. 5 and 8. The recovered shape, shown in 
Fig. 6 is (p,q) = (0.95+JJ;^,0.81+o; ( /o) at 1 R e . Our shape is con- 
sistent with the result from Statler (2001) of p > 0.99 and q ~ 0.7. 
Given that his method is completely different, it is reassuring to see 
that we can even reproduce the shape of the confidence intervals (but 
see van den Bosch et al. 2008). The allowed viewing angles cover a 
large range, but prefer strongly inclined (face-on) views, which is 
also consistent with the results from de Lorenzi et al. (2009). 

It is important to notice that at the 3o" level the shape is not con- 
strained well, allowing almost all viewing angles (f}, (p) and a large 
allowed range in shapes (see Fig. 6). A pure oblate axisymmetric 
spheroid is excluded at the 2a level, and this happens because the 
axisymmetric model cannot reproduce the twist in the zero velocity 
curve (A/ 2 > 200). This is shown in Fig. 5. The differences between 
the axisymmetric (third panel from the left in Fig. 5) and triaxial 
(second panel from the left in Fig. 5) model are not very prominent; 
the most visible change can be seen in the (twisting) shape of the 
zero velocity curve. 

Since S06 uses bi-symmetrized kinematics and the axisymme- 
tric models use different intrinsic mass bins it is not possible to 
directly compare the % 2 of those models with ours. To do a direct 
comparison we recreated the original axisymmetric model from S06 
with the triaxial machinery, without bi-symmetrizing, but with the 
kinematic misalignment correction (see §3.2.1). We expected that 
this original axisymmetric model would fit the data better than the 



axisymmetric model - because the latter does not correct for the 
misalignment - but this is not what we found (rightmost panel in 
Fig. 5). The kinematics of the best-fit triaxial axisymmetric model 
are statistically a significantly much better fit, with a difference in 
A% 2 > 900. The differences show up in the twist of the zero velocity 
curve and the 'hexagonal shape' of the velocity dispersion. It seems 
that in this case the twist of this galaxy was overestimated due to the 
inaccurate determination of the photometric or kinematic PA. Both 
are possible because NGC 3379 is very round, has some isophotal 
twist and does not have a strong velocity field. Luckily, the triax- 
ial modeling (as opposed to the axisymmetric modeling) does not 
depend on these measurements, as it only requires that the relative 
orientation between the photometry and kinematics be known. 

Given that the shape can not be constrained accurately our 
choice of the M. might influence the recovered shape. To test if this 
was the case we checked to see if the recovered shape would differ 
if we set M.= 0. The shape recovered was not different. 

The intrinsic orientation of the central gas and dust disc in this 
galaxy is interesting due to its apparent misalignment of -~50 degrees 
with the main body of the galaxy. The only stable configurations 
in a stationary triaxial geometry are in the principal planes. Statler 
(2001) did a thorough analysis and concluded that, if the disc lies 
in a plane, the only option that he could not rule out was a polar 
ring. Lauer et al. (2005) showed that stellar photometry inside 1 
arcsecond has a sudden PA twist of more then 20 degrees, towards 
the gas disc. The gas disc has a size of 4 arcsec and is thus at larger 
radii than this stellar feature, but the two could be connected. Also, 
S06 found evidence that the gas disc might be warped using an 
ad-hoc model, indicating that a simple stable gas disc in a plane 
might actually not be a good description. Our current mass model 
does not include any isophotal twist and therefore does not predict if 
these two features are intrinsically aligned. To do that, a mass model 
with isophotal twist would be needed. As an added benefit such a 
mass model with isophotal twist will lower the amount of possible 
deprojection, which would indirectly help constrain the shape of 
this galaxy. 

In our current modelling without isophotal twist, the allowed 
range in viewing angles is large, and it might thus be possible to 
place the ring in a principal plane. The ring is misaligned 45° from 
the photometric PA, so to place the ring in a principal plane we need 
a misalignment of the PA of 45°. Our allowed models do include 
these extreme misalignments of \|f = 45 at the 3o level (Fig. 6), but 
the other viewing angles are then quite restricted: at \|/ = 45 only 
(40°< v> < 60° , - 10°< (p < 20° ) are within the 3<T contour, essentially 
disallowing the disc in either the x — y or y — z plane. From our 
modelling the polar ring is the only possibility, as the inclination 
of our best-fit models lies below v}<43°, which is exactly what is 
needed for the polar ring according to Statler (2001). 

3.2.3 The black hole in NGC 3379 

Now that we have a handle on the shape, we investigate whether 
the inferred shape affects the recovered M. . We used the six best- 
fitting mass models, while changing the M. and the M/L. The results 
are shown in Fig. 7. In this figure the shape is parameterized as 
(0.95 — p)/2 + q, which is completely arbitrary, but allows us to 
plot two-dimensional contours and show how the % 2 minimum 
is bracketed. The best-fit shape is independent of the chosen M., 
showing that the recovered shape does not depend on the fixed black 
hole mass that was used in the previous section. 

The best fitting M. is (4 ± 1) x 10 8 M and the M/L is 
(3.0 ±0.2) Mq/Lq./. Surprisingly, this M. estimate from the tri- 
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Figure 5. Stellar mean velocity (Top) and velocity dispersion (bottom) of NGC 3379. From left to right: the SAURON observations, best-fitting triaxial model, 
axisymmetric model and an axisymmetric model without the misalignment. The contours show representative isophotes and the straight line is drawn to 
guide the eye to the zero-velocity-curve. The triaxial model reproduces the data best, followed by the oblate model and last the oblate model corrected for the 
misalignment. 
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Figure 6. Recovered shape and viewing angles for the model of NGC 3379 
with kinematic misalignment. The thin and thick contour represent A% 2 =152 
and 456, which represent 1 and 3 o" confidence intervals. Best-fit model is 
nearly as round as observed and there is no strong constraint on the viewing 
angles. Figure layout identical to Fig. 1. See §2 and §3.2.2 



axial model is more than twice as large as (1.4^'q) x 10 8 from 
the edge-on axisymmetric estimate from S06. To show the qual- 
ity of the models we show the OASIS kinematics and models with 
different M. in Fig. 8. For all the different black hole masses, the 
mean velocity field is reproduced extremely well, but the disper- 




- + + + + + + 



-+ + + + + i , + ~ 

0.65 0.70 0.75 0.80 0.85 

Shape of NGC3379 

Figure 7. The M. (top) and the M/L (bottom) confidence levels as function 
of the different triaxial shapes (see text). The small red crosses represent 
the values for which models where computed and the big black crosses 
indicate the best-fit model, respectively. The contours show increasing level 
of confidence. The inner contour indicates a A% 2 of 42.3, further contours 
represent integer steps of 152, which correspond to the 1,2,3 (thick),4,... a 
confidence intervals on the shape. 



sion is only properly reproduced by the 4 x 10 8 M ( 7, M. model. Our 
result is also significantly above the axisymmetric (near) face-on 
model of Gebhardt et al. (2000, 2.0 x 10 8 M ) and S06. Our M. is 
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Figure 9. Orbital structure of NGC 3379. Top plot shows orbital anisotropy 
6r/(5t = y^o 2 ,/^ 2 , + Og) an ^ tne bottom plot shows the orbit type as 
a function of radius. The model is mildly radially anisotropic outside R. 
and strongly anisotropic inside R.. The short axis tubes (blue solid line) 
dominate the galaxy outside R«, while the box orbits (black dotted line) 
become more important inside R. . The long axis tubes (red dashed line) are 
roughly constant at 20%. 



just outside of the scatter of the M.-a relationship, as that predicts 



10 : 



8.13±0.06±0.27 



(1.4_i'2) x 10 Mq, placing our estimate on the 
heavy side of this relationship. S06 showed that the gas disc inside 
cannot be fit by simple Keplerian motion, implying that the gas 
disk is disturbed and may not be a good candidate for probing the 
dynamical BH mass. Nevertheless they constructed an ad hoc and 
non-unique model and estimated (2.0 ±0.1) x 10 8 M©, which is 
lower than our estimate. 

Even with a mildly triaxial shape, the long axis and box orbits 
can contribute a significant fraction (Hunter & de Zeeuw 1992). 
Inside one intrinsic R e our model consists of 70%, 20% and 10% 
short axis, long axis tubes and box orbits respectively. The orbital 
structure (Fig. 9) of the triaxial model reveals that NGC 3379 is 
radially anisotropic inside the sphere of influence of the black hole 
(R.) and at most mildly radially anisotropic outside. This is different 
from S06, which showed NGC 3379 to be isotropic inside the core 
radius. In our model, the box orbits contribute most of the mass 
inside R., and thus the model becomes strongly radially anisotropic 
in the center (Fig. 9). This is very different from an axisymmetric 
model, in which box orbits cannot exist. 

The box orbits in the center could even be the cause of our 
high M.. In the face-on view of these models (28°< ■& < 49°) the 
stars on box orbits in the center have the highest dispersion in the 
direction perpendicular to the viewer and can therefore affect the 
central observed dispersion. This is exactly opposite to a mechanism 
suggested by Gerhard (1988) that essentially makes the black hole 
unnecessary by viewing the galaxy down the x-axis (end-on) - the 
box orbits would then account for the high velocity dispersion in 
the center 2 . 



4 RELIABILITY OF THE BLACK HOLE MASS 
ESTIMATES 

To get to our best fitting models of M32 and NGC 3379, we had to first 
assume a M., find the best-fitting shape and then find the best-fitting 
M., because the alternative - searching the full parameter space - is 
computationally unpractical. We search the shape parameter space 
first because an initial guess for M. can be done using M.-a and we 
know from van den Bosch & van de Ven (2009) that the influence of 
shape on the quality of the fit is much bigger than that of the black 
hole, usually more than a factor ten in A% 2 . However, this procedure 
is not guaranteed to find the global minimum. 

To ensure that we do find the minimum, we marginalized over 
all the best-fitting shapes when searching for the M.. This showed 
that for both galaxies we found that the best-fitting shape was un- 
changed by the improved M. and that even if we would have chosen 
a different M. beforehand we would still have found the same min- 
imum. Surprisingly, this was also true for NGC 3379, of which the 
shape is not constrained well and the M. estimate changes. This is 
evidence that the recovery of the intrinsic shape and M. are indepen- 
dent, which simplifies future modeling. Figs. 2 and 7 also showed 
the reverse: that the M. estimate does not depend significantly on 
the shape and thus that our M. estimate is reliable, even if we do 
not get the intrinsic shape perfectly correct. The addition of a dark 
halo in the models would be the next step. Gebhardt & Thomas 
(2009) showed that adding a dark halo can increase the M. if the 
constant M/L is overestimated. However this is unlikely to happen 
for NGC 3379 as the best-fit dynamical M/L is already very close to 
the M/L derived from the stellar population (Cappellari et al. 2006). 
The amount of dark matter in this galaxy has recently been con- 
strained by Weijmans et al. (2009), who measured stellar absorption 
line kinematics at large radii. By combining both this outer data plus 
the OASIS observations the effect of the dark matter on the black 
hole mass estimate could now be determined. 

The final test would be to model an analytical triaxial test 
galaxy with a realistic density profile and with a central black hole. 
However, we are not aware of the existence of such self-consistent 
models 3 . As an alternative, galaxies generated using N-body simu- 
lations could be used (similar to Thomas et al. 2007). However we 
shall refrain from doing this test now, as we expect that it would not 
change the results for the two nearly oblate axisymmetric galaxies 
discussed here. 



5 DISCUSSION AND CONCLUSIONS 

We have shown two applications of the triaxial orbit super-position 
method to galaxies that had their central black hole mass measured 
with axisymmetric models. The reason for this was two-fold: First, 
we confirmed that we obtain the same M. when using our new triax- 
ial implementation in the oblate axisymmetric limit, confirming that 
our method is reliable. Secondly, we explored what might happen 
when the assumption of axisymmetry is relaxed. 

For the nearby elliptical galaxy M32 we obtained identi- 
cal results to previous axisymmetric modeling finding a M. of 



3 It is possible to construct analytic triaxial galaxy models with a central 
black hole and an f(E) distribution function, but the isodensity surfaces are 
then identical to the equipotentials, so that the model would be non-consistent. 
The value of using such models as test cases is then limited as we do not 
2 This is also in contradiction with our end-on model of M32, which did expect this arrangement of isodensity surfaces and equipotentials to arise in 

need a black hole. real spheroids, which moreover have anisotropic velocity distributions. 
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Figure 8. Comparison of the OASIS central point symmetrized stellar kinematics (left) of NGC 3379 and models with different black hole masses ( 1,4,7) x 10 
Mo on the right. Stellar mean velocity (Top row) ranges from -60...60 kms~'and the dispersion ranges from 190.. .230 kms -1 . The 4 x 10 s M is the best-fit 
model, while the others are not capable of reproducing the observed central dispersion. The twist in the mean velocity field is reproduced perfectly in the 
triaxial model, which is something that would be impossible for a pure axisymmetric model. While the higher moments are fitted they are not shown as their 
contribution is not nearly as important as the first two moments (They are shown in S06). 



(2.4 ± 1.0) x 10 6 M Q . Our best-fitting shape of M32 is very close to 
oblate axisymmetric, excluding strongly triaxial shapes. The view- 
ing angles are not well constrained. 

We also duplicated the models of NGC 3379. After confirming 
the results from S06 in the axisymmetric limit, we expanded our 
search to include triaxiality. While the shape is not strongly con- 
strained, we found that this galaxy is seen almost face-on and can 
best be described with a triaxial model. The best-fitting models 
are very round in the center and become fairly flattened at larger 
radii. These results confirm the claims about the intrinsic shape from 
Capaccioli et al. (1991), Statler (1994) and Krajnovic et al. (2008). 

The black hole in this triaxial model weighs (4± 1) x 1O 8 M0, 
which is two times bigger than the axisymmetric estimate from S06 
and Gebhardt et al. (2000). We speculate that the difference in the 
estimate is due to the combined effect of changing from an edge-on 
to a nearly face-on view and the inclusion of the strongly radial box 
orbits in the modeling. 

The significance of the change of the M. of NGC 3379 is un- 
clear. While this individual result is still within the scatter of the 
M.-o relation, it can greatly affect all empirical relations based on 
M., like M.-a, if similar effects are seen in more galaxies. If many 
intrinsically triaxial galaxies would have their black hole mass un- 
derestimated using axisymmetric modeling, the empirical relations 
would be systematically underestimating black hole masses at the 
top end, as we believe that triaxial galaxies dominate at the high 
mass end (e.g. Kormendy & Bender 1996; Emsellem et al. 2007). 
To study this in more detail, it is necessary to study the nuclei of the 
most massive galaxies with triaxial models. 

Also, the axisymmetric models and their M. estimates have 
only been tested with purely axisymmetric (and spherical) test mod- 
els. This means that if axisymmetric models cannot accurately re- 
cover black hole masses in triaxial galaxies, and if most galaxies 



are significantly triaxial, then this would also have an impact on 
black hole demography 4 . A strong hint in this direction is given by 
Thomas et al. (2007), who showed that the axisymmetric models 
have difficulty estimating the M/L in triaxial galaxies. Based on their 
results and our NGC 3379 model, we speculate that the M. estimates 
in triaxial galaxies derived using axisymmetric modeling will have 
systematic errors, because the recovery of the M. is strongly linked 
to the M/L. 
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(1) Original (2) First Interpolation (3) Self Calibrated (4) Sys. Correction (5) Spaxels Symmetrized (6) Point Symmetrized 




-15 -10 -5 5 10 1515 -10 -5 5 10 1515 -10 -5 5 10 1515 -10 -5 5 10 1515 -10 -5 5 10 1515 -10 -5 5 10 15 



Figure Al. Example of the point symmetrization process on the SAURON observations of NGC 3377 from Emsellem et al (2004). The top panels show the 
mean stellar velocity (-1 10...+1 10 km s - ') and the bottom row shows the error (2. ..20 kms~'). (1) The original observations. (2) The errors after being expanded 
onto the spaxels and the kinematics after the first interpolation step. (3) After the self calibrating interpolation steps. (4) After the velocity field is corrected for 
the systematics offsets. (5) After the spaxels are point symmetrized. (6) After the final step of binning together the spaxels back into bins. Notice how the errors 
in the bottom region have not changed after symmetrization, because those bins do not have a point-symmetric counterpart. 



APPENDIX A: POINT SYMMETRIZING VELOCITY 
MAPS 

Almost all orbit-based modellers usually symmetrize their input 
kinematics (e.g. Gebhardt et al. 2003; Cappellari et al. 2006), as 
their model assumptions enforce bi- or point symmetry anyway. 
There are numerous reasons to do this. Most commonly this is 
done to reduce the noise in the observations, facilitate x-by-eye 
comparison, force the observations to be bi- or point symmetric, or 
because they want to reduce the number of observables. It is also 
a good method to remove systematic effects, including systematic 
offsets in the odd moments. The symmetrization also improves the 
linear Gauss-Hermite moments reconstruction used in our models 
(Rix et al. 1997). 

Symmetrizing velocity maps is a degenerate problem and there 
are an infinite number of solutions and none are perfect. In this 
appendix we describe a novel method to point- symmetrize SAURON 
velocity fields. It is accurate, conserves the amount of spatial in- 
formation and propagates the errors, without making unnecessary 
assumptions. The IDL-script itself is available in the electronic tar- 
ball on arxiv.org. At the very least, this method is useful to determine 
the systematic offsets in the odd velocity moments. The multi-step 
process is depicted in Fig. Al. 

We make three assumptions: First we assume that the Gauss- 
Hermite moments are orthogonal and un-correlated. This assump- 
tion is also enforced by the dynamical model itself and is therefore 
'fair', but not true. Second, we assume that the velocity field varies 
linearly along spatial coordinates, which is generally true to first 
order. Third, we assume that we can safely symmetrize without 
worrying about the PSF 

The SAURON kinematics are typically binned (Cappellari & 
Copin 2003). This means that most kinematic observations (and 
associated error) span several spaxels (i.e. lenslets). But still every 



spaxel has an individual flux and Signal-to-Noise measurement. We 
assign every spaxel an error based on its relative flux 5 within that 
bin and the error of the kinematic observation (Ajy„) of that bin as 
follows. 



. spaxel 



spaxel 



(Al) 



Where is the total flux in that bin and £fP axel j s me fl ux m the 
spaxel i (2 in Fig. Al). This definition conserves the error of the 
original bin when the spaxels are combined back into a bin later. 
The spaxels can be recombined into the bins using the weighted 
mean: 



* .spaxel .spaxel 



spaxel 



(A2) 



yspaxel ^ t j ie kinematics of each individual spaxel, which are yet 
unknown. To estimate them we construct a linear interpolation of 
the velocity field over the individual spaxels, using the bin centers 
as the nodes in the linear interpolation 6 (2 in Fig. Al). To be as 
conservative as possible we need to make sure that after this linear 
interpolation the combined spaxels still reproduce the original kine- 
matics. Using eqn. A2 we compute for each bin what the difference 
is between the original kinematics and the interpolated one. This 
difference £)/,,„ is then used to self-calibrate the linear interpola- 
tion step, by adjusting the velocities at the nodes using £)/,,„ and 



5 It is also possible to adapt S/N instead of flux. 

6 The exact way in which this interpolation is done is not important. We 
use the default IDL interpolater TRIGRID, which does a decent job of 
extrapolating too. 
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recomputing the interpolation and this is repeated as needed. This 
self-calibration typically converges quickly and ensures that the 
velocity map interpolated over the spaxels reproduces the original 
binned velocity map when binned back (3 in Fig. Al). 

In a triaxial stellar system the odd kinematic moments have to 
be anti-correlated across the center. The SAURON reduction pipeline 
centers the central spaxel on the galaxy center. Therefore it is trivial 
to find the point symmetric counterpart of each spaxel. We can now 
use this information to correct for systematic offsets (e.g. recession 
velocity and template mismatch) in the velocity map of the odd mo- 
ments. It is important to do this before symmetrizing, because bins 
that do not have a point symmetric counterpart will otherwise not get 
corrected for this offset. For each spaxel that has a counterpart, we 
compute their weighted mean velocity and combined error. We then 
take those values plus the value of the central spaxel and compute 
the total weighted average, which is the systematic offset. We add 
this systematic offset to the velocities of the individual spaxels (4 in 
Fig. Al). 

Now we do the actual point-symmetrizing of the spaxels. This 
step is almost identical to the previous step, but has some significant 
differences. For each spaxel that has a counterpart, we compute their 
weighted mean velocity and combined error. We then replace the 
values of those spaxels with the weighted mean velocity that we just 
computed and replace their errors with the combined error multiplied 
by s/2. The error is corrected, because we just stored the velocity in 
two spaxels and we do not want to artificially decrease our errors (5 
in Fig. Al). After we have point-symmetrized all the spaxels (once) 
in this way, we combine the spaxels (and their errors) back into 
the original bins. Now the point symmetrization is complete (6 in 
Fig. Al). The end result has exactly the same bins as the original 
observations and has thus conserved the spatial information. 

In principle it is possible to adapt this routine to bi- 
symmetrization too. However, in the bi-symmetric case no unique 
counterpart of the spaxels exist, because the position angle (PA) of 
the (presumed axisymmetric) galaxy does not have to coincide with 
the pixel grid. To extend our algorithm to bi- symmetrization another 
interpolation or supersampling needs to be added. Alternatively, ap- 
pendix C in Krajnovic et al. (2006) describes a way to bi-symmetrize 
velocity fields, which interpolates the kinematics between the bin 
centroids, without self-calibrating, using an un-weighted mean and 
without propagating the errors. Bi-symmetrization also critically 
depends on the existence and determination of the global PA. Since 
we model (axisymmetric) galaxies here with our triaxial method, we 
can avoid these issues by applying point-symmeterization instead of 
bi-symmeterization. 

The method described here is not perfect and no method for 
velocity map symmetrizing ever will be. In principle it should be pos- 
sible to relax some of our assumptions and design a better algorithm. 
Instead, it would be much more useful to update the orbit-based 
models, so that they are robust against the systematics in the unsym- 
metrized kinematics. 

This paper has been typeset from a TpXJ KTj3£ file prepared by the 
author. 
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